
\documentclass{article}

\usepackage{fancyvrb}
\usepackage{verbatim}
\usepackage[top=4cm, bottom=4cm, left=3.5cm, right=3.5cm]{geometry}
\usepackage{graphicx}

\newcommand{\HRule}{\rule{0.5\linewidth}{0.2mm}}
\newcommand{\seplib}{\texttt{seplib} }
\newcommand{\textind}[1]{\mbox{\scriptsize{#1}}}
\renewcommand{\vec}{\mathbf}

\newenvironment{listing}[1][]%
{\bigskip \noindent \textit{\large{Listing #1}} \verbatim}%
{\endverbatim  \bigskip}

\title{\seplib tutorial - version 0.0.1}
\author{Jesper S. Hansen}


\begin{document}

\maketitle

\clearpage

\tableofcontents

\clearpage

\section{Introduction}
\seplib is a flexible lightweight molecular dynamics simulation library. 
The main motivation behind \seplib is to provide 
scientists with a package that contains the basic routines and tools for 
molecular dynamics simulations 
and which is easily extended at a low level compared to other 
popular packages like GROMACS, NAMD, etc. The library
is inspired by D.C. Rapport's book "The Art of Molecular Dynamics simulations", 
Cambridge University Press, but has more "safe flexibility"
incorporated, as well as different  
ways of initialising the simulation, simulating molecules, different
integrators and so forth.   
 
This tutorial is a part of the library. The reader is referred to the
example programs that  
come with the \seplib package. These are found in the sub-directory
\texttt{prgs/} under the \seplib home  
directory. In appendix \ref{app1} there is a list shortly descriping
each program. 

\seplib features
\begin{enumerate}
\item Fully optimized algorithms for the force calculations
  (neighbour-list and cell-list).  
\item Possibility for the user to implement any pair interaction potentials.
\item Molecular force field including bond-stretching, bending and
  torsion potentials. 
\item FENE bond interactions.
\item NVE, NVT and NPT simulations.
\item Virtual lattice interactions.
\item Simple chemical reactions.
\item and more...
\end{enumerate}

\seplib is released under GNU Public License (GPL) version 3 and is
written using ISO-C99,  
and it is therefore assumed that the user has experience with the C
programming language. Also, to benefit from this tutorial the reader
must have a basic knowledge of molecular dynamics.  

\section{Installation}
To install \seplib you simply use (from the any shell):
\begin{description}
\item[] \texttt{\$ tar zxvf seplib-X.X.X.tgz}
\item[] \texttt{\$ cd seplib}
\item[] \texttt{\$ ./configure}
\item[] \texttt{\$ make}
\item[] \texttt{\$ make test} (optional)
\item[] \texttt{\$ make tools} (optional)
\item[] \texttt{\$ make install} (as root - optional)
\item[] \texttt{\$ make doc} (optional)
\end{description}
\texttt{X.X.X} refers for the version number. 
If the \texttt{make test} procedure was successful 
you can run a simple test program by typing
\begin{description}
\item[]\texttt{\$ ./prg0}
\end{description}
which should run the program \texttt{prg0} and produce some output to
the screen. 
The source for \texttt{prg0} can be found in the \texttt{prgs/prg0.c} source 
file. 

\section{First program}
The directory \texttt{prgs/} contains a set of programs showing how 
simple atomic and molecular systems can be simulated 
using \seplib. Below is a slightly modified version of \texttt{prg0.c}, 
where we simulate a simple 12-6 Lennard-Jones 
system at the initial state point $(T, \rho) = (0.728, 0.844)$ which
corresponds  
to the triple point. Note, we will use reduced molecular dynamics units 
throughout the tutorial.

\begin{listing}[1]
1: #include "sep.h"
2: 
3:  int main(void){
4:    sepatom *atoms;
5:    sepsys sys;
6:    sepret ret;
7:    double t, dt, dens, lbox, etot, sump;
8:    int natoms, n;
9:   
10    // Setting parameter values 
11:   dens = 0.844; 
12:   natoms = 216;
13:   dt = 0.005;
15:   lbox = pow(natoms/dens, 1./3.);
16: 
17:   // Allocating memory 
18:   atoms = sep_init(natoms, SEP_NO_NEIGHB);
19:   
20:   // Setting up the system
21:   sys = sep_sys_setup(lbox, lbox, lbox, 2.5, dt, natoms, SEP_NO_NEIGHB);
22:
23:   // Initializing the positions and momenta
24:   sep_set_lattice(atoms, sys);
25:   sep_set_vel(atoms, 0.728, sys);
26: 
27:   // Main loop 
28:   t=0.0; n = 0;
29:   while ( t <= 100.0 ){
30:	 
31:     // Reset  
32:     sep_reset_retval(&ret);
33:     sep_reset_force(atoms, &sys);
34:
35:     // Evaluate forces acting on between part. 
36:     sep_force_pairs(atoms, "AA", 2.5, sep_lj_shift, &sys, &ret, n, SEP_ALL);
37:
38:     // Integrate particles forward in time 
39:     sep_leapfrog(atoms, sys, &ret);
40:        		    
41:     t += dt; n++;
42:  
43:   }
44:   
45:   // Freeing memory 
46:   sep_close(atoms, natoms);
47:   sep_free_sys(&sys);
48:
49:   return 0;
50: } 
\end{listing}
Lines 4-6 declare three different structure variables. Briefly, the
\texttt{sepatom} type  
contains information about the individual atom (or particle) like the
postion and momentum,  
the type, mass, charge and so forth. \texttt{sepsys} keeps information about 
the system like number of particles and the simulation 
box. The \texttt{sepret} type stores useful data obtained from the
force evaluation and the integrator which enables efficient run-time
data analysis.  

Lines 7-8 declares a few simple helpful variables for later use.

In lines 11-15 the system density, number of atoms, integrator time
step and the simulation box  
lengths are specified; in this particular example we simulate 216
atoms in a cubic box with  
density 0.844. 
  
Line 18 is very important: the function \texttt{sep\_init()} allocates
the necessary  
memory to store the information about the atoms and sets default
values. The default atom mass is   
1 (in reduced units), the default atom type is \texttt{A} and the
default atom charge is 0. The  
second argument to \texttt{sep\_init()} is an integer indicating the
maximum number of particles  
that can be stored in the particle's neighbour-list. 
The neighbour-list is used for large systems and is therefore 
not applied here which is indicated with the macro
\texttt{SEP\_NO\_NEIGHB}. If you want to use the 
neighbour-list you can use the macro \texttt{SEP\_NUM\_NEIGHB} or an
integer value  
usually smaller than 1000. See Appendix \ref{appendix:macros} for macros 
you can use with \texttt{sep\_init()}. 

The system information is stored in line 21: it should be clear what
the arguments  
are, except for argument 4 and 7. The former is the maximum cut-off
length used for  
the interaction potentials; here we only have a single interaction namely the 
Lennard-Jones interaction, so 2.5 is the cut-off for this. The 7'th
argument is an integer  
specifying how frequent the neighbour-list should be updated. Again, we do 
not use neighbour-lists so we can simply use the macro
SEP\_NO\_NEIGHB. If you use the  
neighbour-list values between 5-10 are recommended.

Lines 24 and 25 initialize the atom positions and velocities; here the atoms are 
arranged on a simple lattice and the velocities corresponds to a temperature of 0.728 in 
reduced units. 

After the initial positions and velocities are specified we are ready to integrate the 
equations of motion forward in time which is done in the main loop,
lines 30-42. Lines 32 and 33  
reset the \texttt{ret} variable and the atom force vectors. The
\texttt{sep\_force\_pairs()}  
function evaluates the pair forces between \texttt{A} particles (argument 2). 
Recall, the default type is \texttt{A}. 
The actual interaction is defined through the cutoff, and argument 4, 
the interaction potential.  In this case the cutoff is set to the standard 2.5 and  
the atoms interact via a truncated and shifted 12-6 Lennard-Jones potential.
\emph{Note, the cutoff of this interaction cannot be larger than the cutoff specified
in the function \texttt{sep\_sys\_setup()}}. We will later discuss how
you can implement any  
type of interaction potential yourself. The last argument to
\texttt{sep\_force\_pairs()} specifies  
if interactions need to be excluded, e.g. if two atoms are bonded we
often exclude their van der Waals  
interaction. In this particular case we have non-bonded Lennard-Jones particles, 
so we can use \texttt{SEP\_ALL} indicating that all interactions are
included. See Appendix  
\ref{appendix:macros} for valid macros for \texttt{sep\_force\_pairs()} .

In line 39 the equations of motion is integrated forward in time using the leap-frog 
integration scheme. 

After the main loop finishes the simulation is finished. Lines 46 and 47 
clean up, i.e. here allocated memory is released. Finally, the program exists 
returning zero value to the operating system indicating a successful program execution.

As stated above this program is a modified version of \texttt{prg0.c}; here you will 
also see how to extract simple system properties from the simulation like energies and 
pressure. 

\subsection{Loading postions from an xyz-file}
Often you will be interested in loading the atom positions from a file. With \seplib you 
can import the positions with the \texttt{xyz}-file format. Each line in this file type 
is composed of four columns (or entries): one for the atom type (a character) and 
three for the atom $x$, $y$ and $z$ coordinates (decomal numbers). 
The \texttt{xyz}-file format has a two line  preamble: 
the first gives the number of particles, the second is an optional text 
string~\footnote{See
  \texttt{http://http://en.wikipedia.org/wiki/XYZ\_file\_format}}. In
\seplib  
the optional text string specifies the simulation box dimensions. For example, if we 
wish to simulate 216 atoms in a simulation box with dimensions $6.349
\times 6.349 \times 10.0$  
the first couples of lines could be
\begin{verbatim}
216
6.349000 6.349000 10.0000
A 0.00000 0.00000 0.00000
A 0.94503 0.00000 0.00000
A 1.89010 0.00000 0.00000
...
\end{verbatim}
Note, that the particle postions are \emph{between 0 and the length of
  the box}. To load an  
\texttt{xyz}-file use the \texttt{sep\_init\_xyz()} function, such that Listing 1, 
line 18 could be replaced with 
\begin{verbatim}
atoms = sep_init_xyz(lbox, &natoms, "position.xyz", 'v');
\end{verbatim}
and line 24 is also removed. Note, you still need to set the initial velocities. Using 
this way to initialize your simulation you need not to set the number
of atoms and box lengths since this  
information given in the file. However, you will still need to
initialize the atom velocities.  
\seplib comes with a series of tools that can help you setup relative
complex systems through  
\texttt{xyz}-files. If you wish to save your system such that you can
load it back into your program  at a later stage and continue your situlation 
you should use \texttt{sep\_save()} and \texttt{sep\_load()} functions, where the 
atom momenta is also stored.  
		      
You can set the atom masses and charges using the functions \texttt{sep\_set\_mass()} 
and \texttt{sep\_set\_charge()}. Their prototypes are 
\begin{verbatim}
void sep_set_charge(sepatom *ptr, char type, double z, sepsys sys);
void sep_set_mass(sepatom *ptr, char type, double m, sepsys sys);
\end{verbatim}
and it should be clear how to call the functions.

\subsection{Large system - neighbour list}
Rarely, one performs simulations of small number of particles. For large systems 
(typically more than 500 atoms or particles), you can use an optimized
algorithm for the force evaluations which  
is the most time consuming part of a molecular dynamics
simulation. Basically, this algorithm is a two step  
method: a linked-list in conjuction with nearest neighbour-list method. 
The algorithm will only work for a large 
enough simulations box such that \emph{the simulation box lengths must
  be at least three times the  
size of the maximum cutoff in all three directions}. To use the
neighbour-list method  
you simply indicate to the system setup function how frequently the
neighbour-list should be updated.  
For example, to update the neighbour-list every 10 time step you
replace line 21 in Listing 1 with: 
\begin{verbatim}
sys = sep_sys_setup(lbox, lbox, lbox, 2.5, dt, natoms, 10);
\end{verbatim}
where the last argument to the function specifies the update frequency.

\subsection{Thermostat and barostat}
With \seplib you can perform NVE, NVT and "NPT" simulations. Note, however that 
in the latter case the pressure is controlled via a Berendsen barostat, i.e. 
this does not correspond to the NPT-ensemble, in fact, the
corresponding ensemble is unknown.  
The temperature is controlled via a Nos\'{e}--Hoover thermostat. The prototypes are 
\begin{verbatim}    
void sep_nosehoover(sepatoms *ptr, double temp_desired, double *nu, 
                    double Q, sepsys sys);

void sep_berendsen(sepatom *ptr, double press_desired,  
                   double beta, sepret *ret, sepsys *sys);     	

void sep_berendsen_mol(sepatom *ptr, sepmol *mol, double press_desired, 
                        double beta, sepret *ret, sepsys *sys);	     	
\end{verbatim}
The Nos\'{e}--Hoover thermostat (or reservoir) has a dynamical
variable which is passed  
to the function as a pointer via argument 3 and the reservoir mass is
given by the fourth argument.  
The third argument in the Berendsen barostat is a relaxation
parameter; usually a value between 0.01-1.0 can be chosen. Notice, the
Nose-Hoover thermostat is an atomic thermostat and should be fine for
most purposes. The Berendsen barostat on the other hand has both an
atomic version and a molecular version, the latter should be used for
molecular systems, of course. 


\subsection{Writing your own interaction function}
\seplib provides you with the standard truncated and shifted 12-6
Lennard-Jones potential with  
a cutoff of 2.5 or $2^{1/6}$ (also denoted the WCA potential) where $\sigma=1$ and 
$\epsilon=1$. You can easily define any other pair 
interaction writing your own potential function. The prototype for
such function must be 
\begin{verbatim}
double fun(double rij2, char opt);
\end{verbatim} 
where \texttt{rij2} is the square distance between the two atoms that
interact and \texttt{opt}  
can be either \texttt{'f'} or \texttt{'u'}. If the former option is given as second 
argument the function returns the 
scalar force, the latter option and the function returns the potential
energy due to the interaction.  

As a simple example, consider the purely repulsive potential 
\begin{eqnarray}
\label{eq:pot12}
U(r_{ij}) = 4 \epsilon \left(\frac{\sigma}{r_{ij}}\right)^{12},
\ \ \mbox{for} \ \ r_{ij} \leq 2.5 \sigma ,  
\end{eqnarray}
where $r_{ij}$ is the distance between atom $i$ and $j$ with position 
vectors $\vec{r}_i$ and $\vec{r}_j$, repectively. $\epsilon$ and
$\sigma$ are the usual Lennard-Jones  
parameters defining the energy and length scales. If $r_{ij} > 2.5
\sigma$ the two particles do not interact;  
we will not shift the potential in this example. The force is the negative 
of the gradient of the potential, that is,
\begin{eqnarray}
\vec{F} = - \vec{\nabla} U(r_{ij}) = -(\partial U(r_{ij})/\partial x,
\partial U(r_{ij})/\partial y,  
\partial U(r_{ij})/\partial z),
\end{eqnarray} 
where $x$ is the $x$ component of the vector $\vec{r}_i - \vec{r}_j$. 
Applying chain-rule the $x$ component of the force is 
\begin{eqnarray}
\label{eq:force}
F_x = - \frac{\partial U(r_{ij})}{\partial x} = -\frac{\partial
  U(r_{ij})}{\partial r_{ij}}  
\frac{\partial r_{ij}}{\partial x}.
\end{eqnarray} 
Now, $r_{ij}=\sqrt{x^2 + y^2 + z^2}$ so $\partial r_{ij}/\partial x =
x/r_{ij}$. Substituting this into  
Eq. (\ref{eq:force}) and differentiating the potential interaction
function, Eq. (\ref{eq:pot12}), we obtain 
\begin{eqnarray}
\label{eq:force12}
F_x = 48\epsilon \left( \frac{\sigma^{12}}{r_{ij}^{13}} \right) \frac{x}{r_{ij}} = 
48\epsilon \left( \frac{\sigma^{12}}{r_{ij}^{14}} \right) x
\end{eqnarray} 
and so forth for the $y$- and $z$-components. 

In the \seplib force function the inter-atomic distancies are calculated, i.e. 
$x, y$ and $z$, so you simply need to code the scalar force factor 
$48 \epsilon \sigma^{12}/r_{ij}^{14}$. In reduced units we simply 
use $\epsilon=1$ and $\sigma = 1$, which simplifies the coding further.
This leads to the following definition of \texttt{fun()}

\begin{listing}[2]
double fun(double rij2, char opt){
  double retval=0.0, irij2 = 1.0/rij2, irij12 = pow(irij2, 6); 
  
  if ( opt == 'f' )
    retval = 48.0*irij12*irij2;
  else if ( opt == 'u' )
    retval = 4.0*irij12; 
  
  return retval;
}  
\end{listing}
To use this potential line 36 in Listing 1 is replaced with:
\begin{verbatim}
sep_force_pairs(atoms, "AA", 2.5, fun, &sys, &ret, n, SEP_ALL);
\end{verbatim}
In this way you can also specify any interactions between different particle types. 

\subsection{The \texttt{sepatom} type}
As mentioned above, one of the motivations for making \seplib is to provide 
computational scientists and students with a ligthweight molecular
dynamics library   
which is easy to extend in order to handle specific tasks and nonstandard  
data analysis. Line 4 in Listing 1 declares a pointer to the type \texttt{sepatom}. 
This type is simply a structure containing information about 
each atom's position, velocity, type, charge and so forth. The \texttt{sepatom} 
type is declared in \texttt{sepstrct.h} under the \texttt{include/} directory. The 
most important structure fields are
\begin{verbatim}
typedef struct {
  double x[3];          /**< Position */
  double v[3];          /**< Velocity */
  double f[3];          /**< Force */
  double a[3];          /**< Acceleration */  
  double m;             /**< Mass */
  char type;            /**< Particle type */
  double z;             /**< Point charge */
  int molindex;         /**< Molecule index (what molecule does this atom belong to) */  
	...
} sepatom;	
\end{verbatim}
Data analysis, for example, can now be done by accessing relevant information. 
The following function will calculate the total momentum in the $y$ direction 

\begin{listing}[3]
double eval_momentum(sepatom *ptr, int natoms){
	
  double sum = 0.0;
  for ( int n=0; n<natoms; n++ ) sum += ptr[n].v[1]*ptr[n].m;

  return sum;
}
\end{listing}
Naturally, you can access the other structure fields in the same
manner. To call this 
function from your program, say Listing 1, you can use
\begin{verbatim}
printf("%lf\n", eval_momentum(atoms, natoms));
\end{verbatim}
Note, the type \texttt{sepsys} has a field \texttt{npart}, so you
could also use   
\begin{verbatim}
printf("%lf\n", eval_momentum(atoms, sys.npart));
\end{verbatim}

\section{Simulation of molecular systems}
With \seplib you can simulate molecular systems. The intra-molecular 
force field function 
includes bond, angle and torsion interactions. The total potential 
energy due to intra-molecular interactions excluding possible pair 
potentials is given by 
\begin{eqnarray}
\label{eq:intrapot}
U(\vec{r}) = \frac{1}{2}\sum_{\textind{bonds}}k_s^{(i)} (r_{ij}-l_b^{(i)})^2 + 
\frac{1}{2}\sum_{\textind{angles}}k_\theta^{(i)} [\cos(\theta)-\cos(\theta^{(i)})]^2 + 
\sum_{\textind{dihed}} \sum_{n=0}^5 c_n^{(i)} \cos^n(\phi), 
\end{eqnarray}
where $k_s^{(i)}$ is the spring force constant for bond type $i$,
$k_\theta^{(i)}$ the  
angle force constant for angle force type $i$, and $c_n^{(i)}$ the torsion 
coefficients for torsional force type $i$. To include intra-molecular
interactions defined  
above you can call the \seplib functions \texttt{sep\_stretch\_harmonic()}, 
\texttt{sep\_angle\_cossq()}, and \texttt{sep\_torsion\_Ryckaert()}. 
The prototypes for these functions are
\begin{verbatim}
void sep_stretch_harmonic(sepatom *aptr, int type, const double lbond, 
                          const double ks, sepsys *sys, sepret *ret);
                          
void sep_angle_cossq(sepatom *ptr, int type, const double angle0, 
                     const double k, sepsys *sys, sepret *ret);
                     
void sep_torsion_Ryckaert(sepatom *ptr, int type, const double c[6],  
                          sepsys *sys, sepret *ret);
\end{verbatim}
Notice that the second argument to each of these functions specifies a
particular type such that you  
can have several different bonds, angles and torsion interactions in
the same simulation as indicated  
by the superscript $i$ in Eq. (\ref{eq:intrapot}). This enables you to
simulate complex  multicomponent molecular systems. 

\seplib also supports Finite Extensible Nonlinear Elastic (FENE) potential. 
The function prototype for this is
\begin{verbatim}
void sep_FENE(sepatom *aptr, int type, sepsys *sys, sepret *ret);
\end{verbatim}
Examples are found in programs \texttt{prg1.c} and \texttt{prg2.c}
under the \texttt{prgs/}  
directory. 

\subsection{Setting up the system - the \texttt{sep\_sfg} tool}
The way you specify the molecular structure(s) with \seplib is through
(i) the topology file  
(with extension \texttt{.top}) and (ii) the corresponding \texttt{xyz} file. 
In this example we will see how to setup a system of 
butane molecules based on a united-atomic-unit model, where the methyl
(the two $\mbox{CH}_3$ groups)  
and methylene (the two $\mbox{CH}_2$ groups) are represented by a
single particle with same Lennard-Jones  
parameters. If the Lennard-Jones length scale parameter is $1 \sigma =
3.91$ {\AA} the \texttt{.xyz}  
file for a single butane molecule can be written as 
\begin{verbatim}
4
One butane molecule 
C 0.000 0.000 0.000
C 0.320 0.229 0.000
C 0.640 0.000 0.000
C 0.960 0.229 0.000
\end{verbatim}
In this model the butane has flexible bonds as well as angle and
torsion interactions. All this is specified  
in the \texttt{.top} file in individual sections. A section starts
with a section specifier. For bonds this  
is 
\begin{verbatim}
[ bonds ]
\end{verbatim}
Likewise, the angle and dihedral section specifiers are \texttt{[
    angles ]} and \texttt{[ dihedrals ]},  
respectively. A text/comment line follows after the section
specification, for example,  
\begin{verbatim}
; This is a comment line
\end{verbatim}
The comment line must start with a semi-colon, ;. After the comment
line you specify  
each bond, angle or dihedral depending on the section specifier. 
For the bond section, each bond is described row-wise using four
entries: (i) The first entry gives the molecules index - just put zero here, 
(ii) the second and third entries are the indicies of the two particles 
forming the bond and (iii) the forth entry sets the bond type. For the angles, three particle 
indicies must be specified to define the angle and four for dihedrals. In the case of a butane 
molecule the particle indicies are 0,1,2 and 3, thus the topology file is 
\begin{verbatim}
[ bonds ]
; Bond section 
0 0 1 0
0 1 2 0
0 2 3 0

[ angles ]
; Angle section 
0 0 1 2 0
0 1 2 3 0

[ dihedrals ]
; dihedral/torsion section 
0 0 1 2 3 0
\end{verbatim}
Once you have defined a single molecule from the \texttt{.xyz} and
\texttt{.top} files  
you can generate an entire system's \texttt{.xyz} and \texttt{.top} files 
using the program \texttt{sep\_sfg} which is compiled using 
\texttt{make tools} under the \seplib home directory. 
Assuming that this program is in your search path (e.g. your
\texttt{bin} directory)  
and you have saved the \texttt{.top} 
file above as \texttt{butane.top} and the \texttt{.xyz} as
\texttt{butane.xyz} you type 
\begin{description}
\item[] \texttt{\$ sep\_sfg butane.xyz butane.top 216}
\end{description}
This will produce two new files \texttt{start.top} and \texttt{start.xyz} 
specifying all the atom postions (the \texttt{.xyz} file), all the
bonds, angles and  
dihedrals (the \texttt{.top} file) for the 216 molecules. These two files
can be read directly in your \seplib program using the
\texttt{sep\_init\_xyz()} function  
described above and the function \texttt{sep\_read\_topology\_file()}
\begin{verbatim}
void sep_read_topology_file(sepatom *aptr, sepsys *sys, const char *file, 
                            int npart, char opt);
\end{verbatim}
The molecules' center-of-mass are positioned on a low density simple cubic 
lattice, and you will probably need to compress your system. \seplib
provides you with two  
functions to compress your simulation box: one will perform an
isotropic compression (all spatial 
directions are compressed), the other one will compress the box in one
specified direction 
\begin{verbatim}
void sep_compress_box(sepatom *ptr, double rhoD, double fac, sepsys *sys);

void sep_compress_box_dir(sepatom *ptr, double rhoD, double fac, 
                          int dir, sepsys *sys);
\end{verbatim} 
Here \texttt{rhoD} is the desired \emph{atomic number density} and
\texttt{fac} is the compression  
factor which should be less than one for compression, however, for too
small values and the system will explode. In 
\texttt{sep\_compress\_box\_dir()} the argument \texttt{dir} can be
0,1, or 2 for the $x$, $y$ and  
$z$ directions, respectively. See \texttt{prg1.c} for an example. 

\texttt{sep\_sfg} lets you define as many different type of molecules
you want in a simulation. In general,  the program is called using
\begin{description}
\item[] \texttt{\$ sep\_sfg m1.xyz m1.top num1 m2.xyz m2.top num2 ... }
\end{description}
where \texttt{m1.xyz} is the \texttt{.xyz} file for the atoms 
in molecule one, \texttt{m1.top} the topology file for molecule one,
\texttt{num1} the number  
of molecules of this type, and so forth. Notice that if the bond,
angle or diehdral types are different  
this should be specified in the \texttt{.top} files as well. 

\subsection{Auto-generate a topology file}
For simple molecules you can use the program
\texttt{sep\_ctop} to generate a topology file from an
\texttt{.xyz} file. For example, to generate the topology file for a
butane molecule you use the command
\begin{verbatim}
$ sep_ctop butane.xyz 1 0.39 0.1 1 1.9 0.1 1 3.14 0.1
\end{verbatim} 
which will write the file \texttt{butane.top} specifying the bonds,
angle and dihedral angle. \texttt{butane.xyz} is the single molecule
\texttt{.xyz} file, the next three arguments \texttt{1 0.39 0.1}
specify how many different bonds types to look for (here 1), and
then the length of the bond with some tolerance. You can specify any
number of bond types using
\begin{verbatim}
 <nbtypes> <length bond type 0> <tolerance 0> <length bond type 1>
<tolerance 1> ...
\end{verbatim}
The same arguments are used to specify angles and dihedrals. WARNING,
this tool will only work for relatively simple molecules, always check
the topology file generated
\footnote{It has been tested successfully with 1,7-dimethylnaphtalen,
  but it may fail for more complex molecules}!  
    
\section{Need for speed}
Execusion speed is always a matter of interest to anyone doing MD
simulations. The focus of \seplib is not too give you a top-tuned MD
library, but it does offer optimization in form of neighbour-list
and limited openMP support. We have shown how the neigbour-list
can be used. To use openMP you simply add a call to the
\texttt{sep\_set\_omp()} function
\begin{verbatim}
void sep_set_omp(sepsys *sys);
\end{verbatim}
This will then use the number of threads given by the operating
system. Note, \emph{speed-up will only be achieved for more than 2
  threads, e.g. 4 cores CPU or 4 CPU-SMP machines}. This is due to the fact
that Newton's third law cannot be applied inside the force evaluation
function. Also, the openMP support can only be used together with the
neighbour-list method, i.e. for sufficiently large systems.

Figure \ref{fig:benchmark} shows second per time-step as a
function of number of particles for a simple Lennard-Jones. The
program was run on an AMD phenom quad-core and the library and program were
compiled with gcc 4.4.3 using optimization compiler flag \texttt{-O2}
and \texttt{-fopenmp}.   
\begin{figure}[ht]
  \begin{center}
    \includegraphics[scale=0.4]{doc/latex/sep_benchmark_omp.pdf}
    \caption{\label{fig:benchmark}Performance using openMP support  
      (lower data points) and without using openMP support (upper data points). 
      The inset shows the speed-up defined in Eq. (\ref{eq:speedup}). 
      See text for technical details.}	
  \end{center}
\end{figure}
The inset displays the corresponding speed-up which is defined as 
\begin{eqnarray}
\mbox{speed-up} = \frac{\mbox{cpu time serial}}{\mbox{cpu time parallel}}
\end{eqnarray}
Notice, the theoretical speed-up is 2. With other compilers and
machines the result may differ.

\section{Run-time sampling}

\begin{tabular}{ll}
Sampler                      &     Arguments                 \\
\texttt{sacf}/\texttt{msacf} &                               \\
\texttt{vacf}                &                               \\
\texttt{gh}/\texttt{mgh}     &                               \\
\texttt{hprof}               &                               \\        
\end{tabular}



\clearpage

\appendix

\section{Example programs \label{app1}}
\begin{description}
\item[prg0]{A simple Lennard-Jones system}
\item[prg1]{Butane with flexible bonds (GAFF + Ryckaert-Bellemans parameters)}
\item[prg2]{SPC/Fw water}
\item[prg3]{An "NPT" simulation}
\end{description}

\section{Auxiliary programs}

\begin{description}
  
\item[sep\_sfg]Molecular configuration generator. This program allows you to 
  setup the configuration for a molecular system (including molecular mixtures) from 
  \texttt{.xyz} and \texttt{.top} files for the individual molecules. The files 
  \texttt{start.xyz} and \texttt{start.top} are written and can used as start files. Usage:
\begin{verbatim}
			sep_sfg <m1.xyz> <m1.top> <num1> <m2.xyz> <m2.top> <num2> ... 
\end{verbatim}
\begin{description}
\item[\texttt{m1.xyz}:] {\texttt{.xyz} file for the atoms in molecule one (string).}
\item[\texttt{m1.top}:] {the topology file for molecule one (string).}
\item[\texttt{num1}:] {the number of molecules of this type (integer).}
\item[]	{and so forth.}
\end{description}

\item[sep\_bcc] Creates an \texttt{.xyz} file with atomic positions in a body-centered 
  cubic lattice. The file \texttt{start.xyz} is written and can be used as a start file. Usage:
\begin{verbatim}
			sep_bcc <nx> <ny> <nz> <gabx> <gaby> <gabz> <type>}}
\end{verbatim}
\begin{description}
\item[\texttt{nx}]{Number of atoms in the $x$ direction (integer)}
\item[\texttt{ny}]{Number of atoms in the $y$ direction (integer)}
\item[\texttt{nz}]{Number of atoms in the $z$ direction (integer)}
\item[\texttt{gabx}]{Distance between atoms in the $x$ direction (decimal number)}
\item[\texttt{gaby}]{Distance between atoms in the $y$ direction (decimal number)}
\item[\texttt{gabz}]{Distance between atoms in the $z$ direction (decimal number)}
\item[\texttt{type}]{Atom type (character)}
\end{description}		 

\item[sep\_augxyz] Augment the atom positions given in an \texttt{.xyz} file to another \texttt{.xyz} file.
  The augmentation is in the $z$-direction. The file \texttt{start.xyz} is written and 
  can be used as a start file. Usage:
\begin{verbatim}
			sep_augxyz <xyz-file 1> <xyz-file 2> <dist (optional)>
\end{verbatim}	
\begin{description}
\item[\texttt{xyz-file 1}]{The file to wich the second \texttt{.xyz} file is augmented (string).}
\item[\texttt{xyz-file 2}]{The file wich is augmented, i.e. the positions given in this file is translated
  in the $z$-direction (string).}
\item[\texttt{dist}]{Optional argument that specifies an additional distance between the particles- default is 1.0 
  (decimal number).}
\end{description}	
\item[sep\_ctop] Compose topology file from a (single molecule)
  \texttt{.xyz} file. Usage:
\begin{verbatim}
   sep_ctop <xyz-file> <nbtype> <lb0> <b0eps> ...<natypes> <ta0> <a0eps> 
             ... <ndtypes> <da0> <d0eps> ...
\end{verbatim}

\end{description}

\section{Macros \label{appendix:macros}}

\begin{table}[h]
  \begin{center}
    \begin{tabular}{llc}
      \hline \hline
      &  &  \\
      Macro             &   Description                        & Value  \\
      &  &  \\
      \hline
      &  &  \\
      SEP\_NUM\_NEIGHB    & Default maximum number of neighbours & 1000     \\
      &  &  \\
      SEP\_NO\_NEIGHB     & No neighbours                        & 0        \\
      &  &  \\
      \hline \hline 
    \end{tabular}
    \caption{\label{table:macro_sepinit} Valid macros for \texttt{sep\_init()}}
  \end{center}
\end{table} 

\begin{table}[h]
  \begin{center}
    \begin{tabular}{llc}
      \hline \hline
      &  &  \\
      Macro             &   Description                        & Value  \\
      &  &  \\
      \hline
      &  &  \\
      SEP\_ALL           & Include interactions between all particles & 0  \\
      &  &  \\
      SEP\_EXCL\_BONDED   & Exclude pair forces between bonded particles & 1 \\
      &  &  \\
      SEP\_EXCL\_SAME\_MOL& Exclude intramolecular pair forces & 2\\
      &  &  \\
      \hline \hline
    \end{tabular}
    \caption{\label{table:macro_sepprfc} Valid macros for \texttt{sep\_force\_pairs()}.}
  \end{center}
\end{table}

\end{document}

